. AfASb-7P'&77( f > 

NASA-TP-2776 19880015239 

■ I 

Nonlinear Programming 
Extensions to Rational 
Function Approximation 
Methods for Unsteady 
Aerodynamic Forces 

Sherwood H. Tiffany 
and William M. Adams, Jr. 





NASA 

Technical 

Paper 

2776 

July 1 988 




NASA 

Technical 

Paper 

2776 




3 1176 01416 8034 


1988 


Nonlinear Programming 
Extensions to Rational 
Function Approximation 
Methods for Unsteady 
Aerodynamic Forces 


Sherwood H. Tiffany 
and William M. Adams, Jr. 

Langley Research Center 
Hampton, Virginia 


fVJASA 

National Aeronautics 
and Space Administration 


1 p rnp r rv./ r-r 
1 , i uL 

r ' r / 
•' L L 

: 

.JlIN - u ' - • J 


j LAMGLEY RESEARCH CENTER 
LIBRARY NASA 

1 HAMPTON. VIRGINIA 


Scientific and Technical 
Information Division 



Abstract 

This paper deals with approximating unsteady 
generalized aerodynamic forces in the equations of 
motion of a flexible aircraft. Two methods of for- 
mulating these approximations are extended to in- 
clude both the same flexibility in constraining the 
approximations and the same methodology in op- 
timizing nonlinear parameters as another currently 
used “extended least-squares” method. Optimal se- 
lection of “nonlinear parameters” is made in each 
of the three methods by use of the same nonlinear, 
nongradient optimizer. The objective of the non- 
linear optimization is to obtain rational approxi- 
mations to the unsteady aerodynamic forces whose 
state-space realization is lower order than that re- 
quired when the nonlinear terms are not optimized. 
The free “linear parameters” are determined us- 
ing least-squares matrix techniques. Selected linear 
equality constraints are solved explicitly or included 
implicitly, using a Lagrange multiplier formulation. 
State-space mathematical models resulting from the 
different approaches are described, and comparative 
evaluations are presented from application of each of 
the extended methods to a numerical example. The 
results obtained for the example problem show that 
the number of differential equations used to repre- 
sent the unsteady aerodynamic forces in linear time- 
invariant equations of motion can be significantly re- 
duced (up to 67 percent) from the number required 
by a conventional method in which nonlinear terms 
are not optimized. 

1. Introduction 

The equations of motion of a flexible aircraft con- 
tain unsteady generalized aerodynamic force terms 
which when expressed in the Laplace domain are 
transcendental functions. The availability of effi- 
cient linear systems algorithms for aeroservoelastic 
analysis and design has provided strong motivation 
to approximate the unsteady aerodynamic forces as 
rational functions of the Laplace variable (refs. 1-5). 
Such rational function approximations (RFA’s) allow 
the aeroservoelastic equations of motion to be cast 
in a linear time-invariant (LTI) state-space form, al- 
though the size of the state vector increases due to 
the RFA’s. This increased number of states due to 
the RFA’s is referred to as the aerodynamic dimen- 
sion. There is always a trade-off between how well 
the RFA’s approximate the aerodynamic forces and 
how small the aerodynamic dimension can be kept. 

The RFA formulations in the literature (e.g., 
refs. i-ii) have varying capabilities to perform 
such a trade-off. The purposes of this paper are 
(1) to place existing RFA approaches on an equal 


footing with respect to constraint imposition and the 
optimization algorithms employed and (2) to com- 
pare their performance. 

Currently there are three basic formulations used 
in approximating unsteady generalized aerodynamic 
forces for arbitrary motion using rational functions: 

1. Least-squares (LS) — references 2 and 5 

2. Modified matrix-Pade (MMP) — references 3, 

4, 6, and 7 

3. Minimum-state (MS) — reference 8 

For each of these rational function formulations, 
“best” approximations are determined, in a least- 
squares sense, to tabulated unsteady generalized 
aerodynamic force data. 

An extended LS approach (ELS) was developed 
(refs. 9-11) which includes the capability of enforc- 
ing selected equality constraints on the RFA’s and of 
optimizing the denominator coefficients in the ratio- 
nal functions with a nonlinear, nongradient optimizer 
(refs. 12-14). In this paper, the MMP approach 
of reference 6 and the MS approach of reference 8 
are extended in a similar fashion. The extended 
approaches are referred to as “extended modified 
matrix-Pade” (EMMP) and “extended minimum- 
state” (EMS) approaches. The extensions bring 
them to the same level and flexibility in selecting con- 
straints as the ELS approach of reference 9 and allow 
comparative evaluation using the same nongradient 
optimization techniques. Specifically, 

1. The approach of reference 6 is generalized to 
include the same constraint selection as employed in 
reference 9. 

2. The approach of reference 8 is generalized to 
allow more flexibility in the selection of constraints. 

3. Both are altered so as to use the nongradient, 
nonlinear optimizer used in reference 9 in the con- 
strained optimization of the free “nonlinear parame- 
ters” in each rational aerodynamic approximation. 

4. All three use the same overall measure of rel- 
ative error as a cost function for the nonlinear op- 
timization; the necessity of normalizing the original 
aerodynamic data is essentially eliminated. 

The “linear parameters” are determined by solving 
the equations of the least-squares problem with ma- 
trix techniques. Equality constraints are either in- 
cluded by using a Lagrange multiplier formulation 
of the equations or solved explicitly. Inequality con- 
straints are included in the nonlinear optimization 
process. 

The basic problems involved in approximating 
unsteady aerodynamic forces and methods used to 
solve them are presented in this paper. State-space 



mathematical models resulting from the different ap- 
proaches are described. Results obtained from appli- 
cation of the three extended methods to a numerical 
example are presented. Compared are the actual fits 
obtained, frequency responses for the corresponding 
systems of equations, and the total approximation 
error as a function of aerodynamic dimension. In ad- 
dition some of the difficulties incurred in optimizing 
for the approximation parameters are discussed. 

It is shown that all three extended approaches, 
ELS, EMMP, and EMS, demonstrate marked im- 
provement over a baseline constrained, least-squares 
case in which nonlinear terms are not optimized. 
Furthermore, it is shown that both the EMMP and 
the EMS approximations compare favorably with the 
ELS approximations but with significant reductions 
in aerodynamic dimension. The corresponding aero- 
dynamic dimension of the EMMP approach is ap- 
proximately three-fourths that of the ELS approach, 
and the dimension of the newer EMS approach is 
approximately one-third that of ELS. However, the 
example selected had only one more input creating 
the aerodynamic forces (namely, incremental angle of 
attack due to vertical gust) than the number of gener- 
alized coordinates. For a case in which the difference 
between the number of inputs creating aerodynamic 
forces and the number of generalized coordinates is 
even larger (such as would occur for multiple control 
surfaces with irreversible actuators), the difference 
in aerodynamic dimension between the EMMP and 
ELS formulations would not be as significant. The 
basis for this statement is developed in the paper. 

The software implementation of these RFA meth- 
ods was developed to enhance the aeroservoelastic 
analysis capability of the ISAC (Interaction of Struc- 
tures, Aerodynamics, and Controls) system of pro- 
grams (ref. 15) and exists as an independent module 
in that system. Papers which describe other modules 
may be found by employing ISAC as a key word in a 
library search. 

2. Approximating Unsteady Aerodynamic 
Forces 

This section discusses the problems associated 
with using tabular unsteady aerodynamic force data 
in the study of stability and transient response char- 
acteristics and some methods used to solve them. 
Section 2.1 is a review of the equations of mo- 
tion. Section 2.2 discusses the unsteady generalized 
aerodynamic forces and aeroelastic modeling. Sec- 
tion 2.3 discusses the rational function form used 
to approximate the aerodynamic forces, and sec- 
tion 2.4 presents a discussion of the constraints used 


to improve the reliability of the models in critical 
frequency regions. 

2.1. Equations of Motion 

Unsteady aerodynamic effects as well as struc- 
tural dynamics must be considered when modeling 
flexible aircraft. A commonly employed approach to 
formulating the equations of motion for an elastic ve- 
hicle is based on a chosen number of vehicle vibration 
modes and Lagrange energy equations. Only small 
perturbations from a level equilibrium flight condi- 
tion are considered herein. The perturbed aircraft is 
then represented by a linearized system of equations 
expressed in terms of generalized coordinates £ z (£) 
as explained in reference 16. Through separation of 
variables the elastic deformation of the wing z(x , y, t) 
can be represented by the sum of products of a finite 
set of shape functions zi(x,y) called mode shapes, 
and the £ t *(t) as 


z{x, y,t) = J2 i z i( x ’ y) 6(01 (2-1) 

1 = 1 

The equations of motion (EOM) can then be formu- 
lated in the time domain as 


= q[F^t)+F g (t)]+F 6 (t) (2.2) 


where 

q dynamic pressure 

M, K generalized mass and stiffness matrices 

G damping coefficient matrix where the 

structural damping is modeled as viscous 

matrix of aerodynamic forces due to 
aircraft and control surface motions 

Fg vector of aerodynamic forces due to gusts 

F$ vector of generalized forces due to the 

controls 


F 


F C |F, 


, partitioned matrix of 


aerodynamic forces 


Equation (2.1) when transformed into the 
Laplace domain can be written as 


[Ms 2 + Gs + K-gQ e ]£(s) 

= QQg<Xg + FhmThm (2-3) 


2 



where 

Qe matrix of aerodynamic force coeffi- 

cients due to aircraft and control sur- 
face motions, equal to [Qy (s)J, for 
i = 1, and; = 1, 

Qg matrix of aerodynamic forces due 

to gusts, equal to [<2y(s)], for i = 

1, . . . , and j — + 1, . . . , n^ + n g 

(Xg nondimensional gust velocities, 

{wg/u v g /u} T 

F H m x ^8 matrix of modal coefficients 
converting hinge-moment outputs to 
generalized forces 

Thm vector of hinge moments output by actu- 
ator (note that Thm is typically a func- 
tion of actuator characteristics, backup 
structure stiffness, and control law feed- 
back; for this paper, identification of the 
additional dynamics and interconnections 
is suppressed) 

Q = j Q g , matrix of generalized aero- 
dynamic force coefficients 

It is the matrix Q of force coefficients which is of 
interest in this paper. Each element Qij of Q is 
defined as 

Qii (») = // APi(X ' V ' 3) z i (x, y) dS (2.4) 

s 

Here A Pj(x,y,s) is the pressure difference due to 
motion in the jth generalized coordinate. 

Rational function approximations (RFA) allow 
the aeroservoelastic equations of motion to be cast 
in an LTI state-space form for which a large number 
of efficient linear systems algorithms are available: 

X = AX + BU (2.5) 

where 

X vector of states 

U vector of inputs 

A, B matrix multipliers 

Thus there is strong motivation for approximat- 

ing the unsteady generalized aerodynamic forces 
(GAF’s) as rational functions of the Laplace variable 
(refs. 1-5). Subsequent sections describe the steps 
required to obtain the RFA’s. 


2.2. Generalized Aerodynamic Forces and 

Aeroelastic Modeling 

This section sets the stage for obtaining the 
RFA’s of the generalized aerodynamic forces by stat- 
ing how the GAF’s are currently generated, by indi- 
cating the condition for which the GAF’s are desired, 
and by introducing the approach that is employed to 
approximate the desired representation. 

2.2 J. Oscillatory motion. Using such techniques 
as doublet-lattice or kernel function theories (refs. 17 
and 18), programs currently available for genera- 
tion of unsteady GAF’s can compute the forces only 
for purely oscillatory motion over a range of speci- 
fied values of reduced frequency. Consequently, the 
GAF’s are not defined as explicit functions of A;, but 
are defined only as tabular functions and are tran- 
scendental. Using these tabular functions, iterative 
methods must necessarily be used to determine eigen- 
values of the flexible aircraft equations of motion. 
These iterative methods tend to be costly, and the 
solutions are exact only for purely oscillatory motion. 
Most importantly, LTI methods cannot be used. 

2.2.2. Arbitrary motion. In order to obtain so- 
lutions in the Laplace domain for both growing and 
decaying motion (s off the 2 cj-axis), it is necessary 
to express the forces as functions of s for the en- 
tire complex s-plane, or equivalently for the non- 
dimensionalized complex p-plane. In lieu of develop- 
ing new aerodynamic theory, the concept of analytic 
continuation is often used to justify extending these 
functions, which are defined only on the frequency 
axis ( s = iu), to the entire complex plane by finding 
analytic functions which agree with the aerodynamic 
forcing functions at all values of frequency. However, 
there are only a finite number of frequencies at which 
tabular data are available; hence, this process is at 
best an approximate analytic continuation into the 
region near the portion of the axis containing the 
tabular data. Since phenomena such as flutter occur 
for points in the complex s-plane which lie along the 
ito- axis, approximations into the region near the axis 
are sufficient for most studies. 

Figure 1 depicts the approximating process, 
where 

k reduced frequency 

p nondimensionalized Laplace variable 

Qij{ik n ) reduced-frequency domain tabular data 
(identified by the open circles) 

Qij{ik) approximating curve, Qij{p) for p = ik 
(corresponding to the solid line) 
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Qij{ik n ) points along the approximating curve at 
reduced frequencies k n corresponding to 
the tabular data (solid circles) 

Sij(ik n ) approximation error between two corre- 
sponding points (denoted by an arrow 
between points) 

2.3. Rational Function Approximations to 
Generalized Aerodynamic Forces 

The most common form of the approximating 
functions used currently for each generalized force 
coefficient Qij of Q is a rational function of the non- 
dimensional Laplace variable p. Each can be ex- 
pressed in the following partial fraction form, where 
til is the number of partial fractions (referred to 
herein as the “order of fit”), which is equivalent to 
the order of the overall denominator polynomial: 

Qij ( p ) = {Ao)ij + {A\)ijp + (A2)ijP 

n L 

+ <M> 

t— 1 C 

which can be rewritten as 

Qij( s ) = (^0 )ij + (A\)ijS + ( A.2)ijS 

UL - Q 

+ < 2 - 7 > 

where 

c 


be = T b e 

Aq = Aq 


^£+2 — ^£+2 

The partial fractions are commonly called lag terms 
because each represents a transfer function in which 
the output “lags” the input and permits an ap- 
proximation of the time delays inherent in unsteady 


aerodynamics (ref. 5). Because tabular data are 
determined for specified values of reduced frequency 
k n the are actually defined only for values of the 
nondimensionalized Laplace variable p = ik n . Vari- 
ations in the matrix form of the rational functions 
result in three currently used approaches, mentioned 
in the introduction, to approximating the unsteady 
aerodynamic force coefficients. 

2.4. Constraints 

It is often desirable to impose constraints on the 
approximating functions. For instance, one might 
wish that the approximating functions agree with 
the tabular data at the steady-state conditions (e.g., 
Qij — Qij at k = 0). Without constraints imposed, 
as depicted in figure 2(a), the fit at zero frequency can 
be poor, although the overall shape of the approxi- 
mating curve may be good. Figures 2(a) and 2(b) 
depict the “steady-state area” (within a shaded box) 
where constraints might be imposed. In some in- 
stances, agreement of the slopes of the curves at a 
specified frequency might be necessary. For criti- 
cal flutter modes, it might be desirable to impose 
agreement between computed data and the approx- 
imating functions near the flutter frequency. When 
constraints are imposed, there is a corresponding loss 
in degrees of freedom for the least-squares solution. 
The possible change in the approximating curve due 
to this loss is also depicted in figure 2(b). 

The standard method of improving the fits, espe- 
cially when constraints have been imposed, is to in- 
crease the order of fit used in the approximation (i.e., 
make n i in eq. (2.7) larger) or to decrease the fre- 
quency range of the approximations. In recent years, 
methods have been developed to improve these RFA’s 
by optimally selecting certain nonlinear coefficients 
via nonlinear programming techniques (refs. 6-9). 
These studies indicate the advantages of optimizing 
the nonlinear coefficients not only to reduce the error 
for a given order of fit, but also to provide a way of 
reducing the order of fit required to achieve a speci- 
fied level of error. These concepts are explored and 
extended in this paper. 

3. Matrix Formulations of Rational 
Function Approximations and Their 
Corresponding State-Space Equations 

There are several variations on the matrix form 
of the rational function approximations for the un- 
steady aerodynamic force coefficients corresponding 
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to the LS, the MMP, and the MS formulations. Each 
matrix formulation results in different aerodynamic 
state vectors. The mathematical models for each 
of the three main variations on the matrix form of 
the rational functions and their corresponding state- 
space equations are defined in this section. Sec- 
tion 3.1 presents the equations for the least-squares 
RFA formulation as developed by the authors and 
used in the comparisons presented in this paper. Sec- 
tion 3.2 defines the equations for the modified matrix- 
Pade RFA formulation, and section 3.3 defines those 
for the minimum-state RFA formulation. 


3.1. Least-Squares (Column-Independent) 

RFA Formulation 

Historically, Roger (ref. 2) and, later, Abel (ref. 5) 
formulated the rational function approximations so 
as to use the same denominator coefficients for all 
the elements Q{j of the matrix Q in order to reduce 
the number of aerodynamic states in an LTI form 
of the equations of motion. In effect, therefore, the 
formulation was “element-independent” with respect 
to the denominator coefficients. It is, nevertheless, 
referred to here as “column-independent” to distin- 
guish it from the MMP approach of section 3.2. 


The condition that each element Qij of Q has the same constant denominator coefficients, or lag coefficients, 
allows the “per element” expressions of equation (2.6) to be valid for the full matrix Q: 


Q(p) = 


Q tip) | Qg(p) 




= A 0 + A !P + A 2 p 2 + ^i+2—j; 


e=i 


p + k 


which, when expressed in the Laplace s-domain, is: 


(3.1a) 


QM = 


Q*M | Q g(s) 


n L 


— Ao + Ais + A 2 s 2 + T, A^ + 2 — ~t~ 
e=i s + b t 


(3.1b) 


Each matrix in equations (3.1) is rectangular, with dimensions n £ x (n^ + n g ). To develop a state-space 
realization for equation (2.3), Q(s) is approximated by Q (s), as defined in equation (3.1b), giving 


Ms 2 + Gs + K - q 


nL 


( A o)f + (Ai)^s + ( a 2)^s 2 + y^(A£ + 2 )f— 

e=i s + k 


11 


CM 


= F HM T HM + q 


nL 


(Ao )g + (Ai) ff s + {A 2 )gS 2 + nw. 




8 + h 


a 9 


(3.2a) 


If the aerodynamic states are defined by 

= TT~ ( ^ 1 (For £ = 1, . . . , ni) (3.2b) 

8 + be ( Oig J 

there would be (n^ + n g ) x ni aerodynamic states in the LTI state-space form (eq. (2.5)). 

3.1.1. LS RFA aerodynamic dimension reduction. An alternative formulation of the aerodynamic states 
described below significantly reduces the number of aerodynamic states when n g is nonzero. This modification 
uses a definition for the aerodynamic states which combines the elastic and gust modes. 
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Equation (3.2) is equivalent to 


n L 


= 1 
2 


{[M - 9(A 2 ) e ] 5 2 + [G - qiAJtls + [K - q( A 0 ) $ ]} 1(a) - q £ 

= + Q [(Ao) ff + (Ai) ff s + (A 2 )gS 

Hence, if the augmenting aerodynamic state vector X a is defined by 

x a = [x ai x a2 ... x a „ L ] : 

where each 


( A £ +2 )£ I (A e+2 )g 


— * 
S + bp I &g 


(3.3a) 


X af = 


(Am- 2 )$ | (A^ +2 ) s 

I 


s + 1 


W- 


= 7Tw (A m) ” = ttt/ 1 


and 


’■*U} 

Ve = (A£ +2 ) t} (For l = 1, ... , n L ) 
then an LTI state-space realization for equations (3.3) is 

sX 0 , = -6^IX a , + I r)' e 


(3.3b) 


(3.4) 


or more explicitly 


sX a , — -b(IX. ae + (Af +2 )^ si + (A£ +2 ) ff sa g 


(3.5) 


Since A^ +2 is of size x (n^ + n g ), each r)' e is of size n £, so there are exactly states for each i — 1, . . . , n^. 
This implies that the number of added aerodynamic states for state-space realization of equation (2.3) due to 
the partial fractions (or lag terms) in equations (3.1) for this column- independent formulation is 


n a = n^riL 


(3.6) 


The state equations for the aeroelastic system for Roger’s RFA formulation are obtained by substituting 
equation (3.3b) into equation (3.3a), converting to first order form, and then combining with equation (3.5). 
The resulting set of LTI state-space equations are 
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5 < 


e 

si 

{ x a 


> = 


T —1 


I ! 


! 0 


0 I M-g(A 2 ) e I 0 


-(K- 9 (Ao) € ) 


+ 


1 

0 

0 

°_! 

M-q{A 2 )z 

0 

1 

0 i 

0 

I 


-1-1 


“(G - q{Ai)^) 


?[I 


0 

1 

1 

1 

1 

1 

( A 3)$ 

-M 

0 

1 

1 

1 

(^n L +2)f 

i 

j o ... - 

1 

0 j 

I 

i 

! o 

I 0 

f hm j 

g(Ao )g 

| 9( A i) ff 

| 9( A 2)s 

1 

0 j 

0 

! ( A 3 )g 

j 0 

1 

0 J 

0 

(A( ni -+2))g 

1 0 


si 

l Xq J 




f T HM 1 


SOfg 

2 

S a„ 


(3.7) 


This formulation was used for most of the analyses reported in references 19-21. 


3.1.2. LS RFAfor irreversible actuators. For the column-independent LS formulation, a significant reduction 
in the aerodynamic dimension is possible if the actuators driving the control surfaces can be considered 
irreversible. In this case, the n s i inputs associated with control deflections are approximated as not being 
affected by aerodynamic and inertial cross-coupling hinge moments. Thus 


6' = T A (s)6' c (3.8) 

my x 1 vector of commanded control inputs 

diagonal transfer matrix relating control deflections to commanded control 
inputs 

Under this approximation, only the first n^, equations of (2.3) are needed since 

(39) 

and hence the rows defining the n^t deflections can be deleted, and terms in the first n^i equations involving 
6' can be expressed in terms of commanded inputs. Thus the aerodynamic dimension for this case, using the 
LS formulation, is only 

n a = n^ni, (3.10) 


where 

% 

T A (s) 
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3.2. Modified Matrix-Pade (Column-Dependent) RFA Formulation 

In the MMP approach, the denominator coefficients b( are fixed for all elements of a particular column, but 
their number and values may vary between columns; that is, they may be “column-dependent.” Thus, for the 
MMP formulation, the per element expression (eq. (2.6)) holds only per column, not for the entire matrix Q 
as in equation (3.1): 




Q j(p) = ( A o)y + (Ai)yp + (A 2 ) y p 2 + E(Am), 


i=l 


p + 


The Q j(p) are then combined as 


Q(p) = Aq + Aip + A 2 p 2 + [Wi(p) 


Mp) 


W n^+n 9 (p)] 


0 


0 


P 



(3.11) 


(3.12a) 


where 


L j L j 

w,-(p)=2(A< + 2)y n (p+m 

«=i k = l 

k±l 


(3.12b) 


and 


n L, 


h j(p) = YliP + hj) 
£=1 


(For j = 1, 


n S 


n„ 


(3.12c) 


If the augmenting aerodynamic state vector, 


is defined by 


X a = [X ai x a2 ... x an$+n / 



' 0 

1 •• 

0 ' 


/°\ 

«X aj = 

0 

0 •• 

1 

x a .+ 

0 


—hji 

0 •• 

~ h jn L - _ 


u J 


= H jX aj + V jT) (For j — 1, . . . , ri£ + ng) 


(3.13a) 


(3.13b) 
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then the state-space equations for this column-dependent formulation are given by 




r c ) 

*1 

l Xa J 


> = 


i ! 


0 ! M-q(A 2 )i ! 0 


0 ! 



1 

0 

l 

/ 

0 

1 

! 0 

-l 

-[K- 9 (Ao) C ] | 

— [G — qr(Ai)^] j 

9 [Wi ... W n i+ n g ) 

1 





1 

! 0 


0 

Vi ••• 0 j 

Hi ... 0 0 

1 

! r - 


: i 

: i 

o ... Vn« 

0 ... H n? 0 


l 

0 1 

1 

0 ... 0 J 

0 ... 0 H n £ +riff ^ 


r J i 

IX 0 ) 


+ 





-l 

0 

I 

0 

0 


0 

L_° 


Fhm 

<?(Ao)0 


j l{M)g 


M-9(A 2 ) € 



0 

0 

0 

0 

I 

° 

1 

! 0 

1 

! 









0 

0 

V 3 

i 0 


' T H m ' 

a g 

< ► 

sag 

< s 2 a g j 


( 3 - 14 ) 


where, again, £ is the vector of Laplace transformed generalized coordinates, and the bar over the Wy and Hy 
denotes the corresponding matrices in terms of s; specifically, 


[W(«)]j 


pftWly 



Since each H ? is a square matrix of dimension n j J . x ., the number of aerodynamic states for this column- 

J 3 3 

dependent MMP formulation is 


n^Jlg 

n a = n Lj ( 3 - 15 ) 

3=1 

If each rii . = n^, then the number of augmenting aerodynamic states for the MMP formulation is rigni 
larger than the column-independent LS formulation. If n ^ actuators can.be considered irreversible, the number 
of augmenting aerodynamic states for the MMP formulation is larger than the LS by an amount (n^/ J rn g )ni l . 
This amount can be sizable. In every case, however, the MMP fits should be as good as or better than the LS 
for a fixed n l because the denominator coefficients are optimized on a per column basis in the MMP approach. 
The power of the MMP formulation for reduction in added aerodynamic states is maximized by allowing both 
nj^. and {6^.} to vary per column. 
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If the lag coefficients are the same (same values for { b ^ } and n^. — rti) for all columns in the MMP 
formulation, then the column-dependent state-space equations can be modified to have the same aerodynamic 
dimension as the column-independent LS formulation as follows. Since each h j(p) = h(p), equation (3.12a) 
can be rewritten as 


r m ••• 0 


Q(p) = Aq + Aip + A 2 p 2 + [Wi(p) ... W n?+ n 3 (p)] 


p (3.16a) 


hfp) J 


= A 0 + Aip + A 2 p 2 + ^j[Wi(p) ... W n? +n s (p)]P 


(3.16b) 


where h(p) is of order n^, and W(p) is of dimension x (n^ + n g ). 
Since [l//i(p)]W(p)p can be expressed in partial fraction form as 


1 

h{p) 


n L 

W(p)p= Yj A M2 


e=i 


p 

p + h 


equation (3.16b) is equivalent to equation (3.1a); hence, the aerodynamic dimension corresponding to the 
state-space realization can be reduced as in section 3.1.1 to that of equation (3.6): 


n a — n^riL 


(3.6) 


3.3. Minimum-State RFA Formulation 

Karpel suggested in reference 8 that in order to find a minimal augmented state vector, the process of 
finding aerodynamic approximations and then determining the state-space equations be reversed. 

He suggested starting with the general form of the state-space equations defined by 
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where 


Mi ••• 0 1 


R = 


L 0 ... — 6jvJ 


(3.18) 


is a diagonal matrix of aerodynamic roots b( of some order N necessary to accurately represent the aerodynamic 
forces. This implies that the aerodynamic states can then be expressed as 


sX a 



(3.19) 


which can be rewritten as 


x Q = [si - it] 1 



(3.20) 


Since R is of dimension N x N, the aerodynamic dimension for this minimum-state formulation is 


n a = N 


(3.21) 


This aerodynamic dimension N is a parameter which must be determined. Karpel (ref. 8) showed that 
finding a minimal augmenting state vector satisfying equations (3.17)-(3.21) is equivalent to finding matrices 
Aq, Ai, A 2 , R, D, and E so that the aerodynamic approximations have the form 


Q(s) = Ao + Ai s 4- A 2 S 2 + D[sl - R] X E s 
which can be written in terms of p as 


so that 


Q(P) = A 0 + Ai 




I - R 



(3.22) 


(3.23) 


Q(p) = Aq + Aip + A 2 P 2 + D[pl - R] X E p (3.24) 

where D is of size n^xN, R is a diagonal matrix of size N, and E is of size N x (n^+n g ). An iteration method 
is required to determine a converged solution for the matrices D, E, Aq, Ai, and A 2 in equation (3.24). A 
higher level nonlinear optimization is designed to optimally select the R. 


3.4. Comparison of Aerodynamic Dimension 

Typically (but not necessarily) , for a given level of total approximation error, 


max{n^.} < JV_ 


n a 

for MS 


n^+Tlg 

< n L .< 

3 = 1 


n a 

for MMP 


n^n L 

n a 

for LS 


(3.25) 


The character of the matrix corresponding to the lag terms and 
required for state-space realization are summarized as follows: 


the resulting aerodynamic dimension 
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Aerodynamic 

Method 


Character of lag terms 

dimension 

Extended 

Common lag coefficients 

in denominator of 

ngn L 

least-squares 

each matrix element 






n^+rig 

Modified 

Different number of and values for lag 

7) n L 

matrix-Pade 

coefficients in each column of Q 

j= i 

Minimum-state 

Elements of 


N 



'ik • 0 




D 

0 • • • s+b N J 

E 



Thus the MMP formulation gives more flexibility than the LS approach in the choice and number of lag 
coefficients and allows them to be determined on a per column basis. In the MS approach, the fact that the 
number of additional states is independent of n £ provides the potential for a large reduction in aerodynamic 
dimension. 


4. Extended RFA Formulations 

In most of the applications of aerodynamic ap- 
proximations in the literature, the denominator coef- 
ficients b£ in equation (3.1) have been specified a pri- 
ori over the range of desired frequencies based on 
engineering judgment. An alternative approach is to 
optimize these nonlinear parameters to improve the 
approximations by improving their selection. Opti- 
mization of these parameters has been included in the 
studies by Dunn (refs. 6 and 7), Tiffany and Adams 
(ref. 9), and Karpel (ref. 8). 

In reference 6, the modified matrix-Pade approx- 
imation formulation (which differs slightly in form 
from the MMP formulation included herein) is in- 
directly constrained to nearly match certain data 
points by applying frequency-dependent weighting 
factors to the corresponding terms in the least- 
squares solution for the linear coefficients. The non- 
linear optimization method of feasible directions of 
reference 22 is employed therein, using an analytical 
scheme to compute gradients, when determining the 
nonlinear parameters. 

In reference 8, employing the minimum-state 
formulation, three fixed equality constraints are 
imposed on all elements of Q. Explicit inclusion of 
these constraints reduces the number of degrees of 
freedom in the least-squares solution. A modifica- 


tion of Davidon’s variable metric method (refs. 23- 
24) and finite differences to approximate gradients 
(ref. 25) are used in determining the elements b £ of 
matrix R in equation (3.24). 

In reference 9, the authors extended the least- 
squares approach (referred to as ELS) to include 
some capabilities which can have advantages over the 
modified matrix-Pade approach of reference 6 and 
the minimum-state approach of reference 8. First, 
the flexibility in selection of the equality constraints 
of reference 11 was included in order to improve the 
fit to the tabular aerodynamic force data in critical 
regions by imposing certain constraints, or to im- 
prove the fit overall by relaxing certain constraints. 
Second, it employed a numerically stable, non- 
gradient, nonlinear optimizer to determine the b i in 
equation (3.1). This algorithm avoids some numer- 
ical difficulties which may arise in the optimization 
process due to numerical computation of gradients 
and due to nearly linear dependency of active con- 
straints. Side constraints on the b £ could also be 
imposed to ensure system stability. 

In order to comparatively evaluate the three RFA 
matrix formulations, the MMP and MS formulations 
are extended in this paper in a fashion similar to the 
ELS formulation. These extensions include the same 
flexibility in selecting equality and side constraints 
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and include the use of the same nonlinear optimizer. 
The extended formulations are referred to as the 
extended modified matrix-Pade (EMMP) and the 
extended minimum-state (EMS) methods. 

In this section, the selectable constraint options 
will be defined; a common measure of performance 
of the RFA’s will be specified; and the algorithms for 
determining the optimum values of the free parame- 
ters will be specified. 

For a given aerodynamic dimension, the opti- 
mization of the nonlinear parameters for the ELS, 
the EMMP, and the EMS formulations results in a 
smaller total approximation error for suitably con- 
strained approximations than for the LS formulation. 
Alternatively, optimization of the nonlinear param- 
eters can reduce the aerodynamic dimension for a 
specified level of total approximation error. For ELS 
and EMMP, the reduction is indirect since the aero- 
dynamic dimension is a function of both and ei- 
ther or + n g , respectively. The EMS form, 
however, is structured in such a way that the aero- 
dynamic dimension equals the order of fit, thus allow- 
ing for minimization of the aerodynamic dimension 
directly. 

4.1. Unconstrained Linear Optimization 

To determine the linear parameters in the RFA’s, 
a measure of error between the approximating curve 
and the actual tabular data for each aerodynamic 
force element is used as an objective or cost func- 
tion. This particular objective function essentially 
eliminates the necessity of normalizing the original 
aerodynamic data. While not apparent in the linear 
optimization, normalization is essential in the non- 
linear optimization (described later) in order for the 
overall objective function, which is the sum of ele- 
ment errors (e^), to be insensitive to differences in 
magnitudes of the aerodynamic force coefficients. 

A square error function is defined by 


En 1 Qijiikn) ~ Qij{ik n ) | 2 

(4.1a) 

^ - Mij 

where 

M ij = max{l, | Qij{ik n )\ 2 } 

(4.1b) 


and {fc n } is a set of reduced frequencies for which 
tabular data are available. 

Each term in the sum in equation (4.1a) is a 
measure of relative error if the maximum magnitude 
of Qij(ik n ) is greater than 1, but is an absolute 
error for magnitudes less than 1. This error function 
essentially normalizes the aerodynamic data prior to 
the nonlinear optimization. 


A rational approximation of the form of equa- 
tion (2.6) is linear with respect to the coefficients 
( A m )ij . If the coefficients (A m )ij are used as de- 
sign variables for a particular element of Q and equa- 
tion (4.1a) satisfies the minimum condition that all 
partials with respect to each design variable be 0, 
that is, 

ds 

— = 0 (For m = 1, . . . , n L ) (4.2) 

then the resulting system of equations is linear and an 
algebraic solution for the coefficients (A m ) 2 y is pos- 
sible. Henceforth, this will be referred to as a linear 
optimization. Even though the objective function is 
quadratic, linear algebraic methods may be employed 
to determine the optimum solution exactly. 

4.2. Constrained Linear Optimization 

Several types of equality constraints are consid- 
ered in references 9-11. Those which have been in- 
cluded here are listed in the next section. 

4.2.1 . Selectable linear equality constraints. 

The types of linear equality constraints from which 
selection can be made (refer to appendix A in ref. 11) 
are 

1. To constrain the values of the RFA’s at 
zero frequency to be the same as the 
values of the tabular data. 

2. To constrain the slopes of the approxi- 
mating functions at k = 0 to be specified 
values or to satisfy a specified relationship 

3. To null out specific coefficients in the 
approximations 

4. To constrain the values of the RFA’s at 
specified nonzero frequency, for which 
tabular data are interpolated if necessary. 

Each constraint is imposed on the entire column j of 
Q and can be expressed in terms of linear equations 
involving the design variables ( A m )ij for a specified 
j as 

C]{Am}ij = C° (For i = 1, . . . , n { ) (4.3) 


4.2.2. Lagrange formulation. Equality constraints 
of types 1 and 3 are included explicitly. The Lagrange 
multiplier formulation for including constraints of 
types 2 and 4 is employed to redefine the fit error for 
a given element of the aerodynamic force coefficient 
matrix as 
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_ En \Qij{ ik n) ~ Qjj{ik n )\ 2 
+ 2Xf j [cJ{A m } ij -Cl] (4.4a) 
= Sij + 2A [Cj - Cl] (4.4b) 


in appendix A. This nongradient optimizer has been 
found to be numerically stable and to possess good 
convergence properties in the applications to which 
it has been applied (refs. 9 and 21). 

The nonlinear parameters are selected via the 
nonlinear optimization algorithm to reduce the total 
column approximation error, 


where Mjj is defined by equation (4.1b). 

The number of constraint equations n c j varies per 
column. The corresponding gradient conditions 


jj = 



Wi 

13 


(For EMMP) 


(4.7a) 


de c - 

——^— = 0 (Form = 0, ..., n L +2 {orn L +2)) 
uyslmjij J 

(4.5a) 

and 

de c -. 

a/x \ = 0 (For m = 1, . . . , n e ) (4.5b) 

define a system of linear equations whose alge- 
braic solution determines the numerator coefficients 
( Am)ij in equations (3.1). (3.11), and (3.24) and is a 
least-squares error approximation which exactly sat- 
isfies the imposed constraints. The D ^ and E^. in 
equation (3.24) are also determined in this process 
since 

Ee = &) » } (For m = i + 2; l = 1, . . . , n L ) 

(4.6) 

for alternate steps in the iteration process converging 
to a solution for D and E (refer to section 4.3.3). 

4.3. Multilevel Optimization: Linear and 
Nonlinear 

Improvements in the approximations can be made 
by increasing the number of lag terms. As shown in 
section 3, however, this increase in the order of fit ad- 
versely affects aerodynamic dimension and increases 
the number of equations required to define the aero- 
elastic system. The fit can also be improved by re- 
ducing the frequency range over which the fits are 
required, but this reduces the range of applicability 
(or robustness) of the approximation. The additional 
free parameters in all the approximations, namely the 
b(> can be optimized to improve the approximations. 
Historically, these parameters, or their equivalent in 
the matrix-Pade approximants, have not been in- 
cluded as design variables since the resulting gradient 
equations, unlike equations (4.5), are not linear. The 
technique employed herein to optimize these non- 
linear parameters is a nongradient method described 


summed over all rows, or the total matrix approxi- 
mation error, 



(For ELS and EMS) (4.7b) 


summed over all columns. The weighting factors W{j 
are used to force some of the elements to have more 
priority than others in determining the nonlinear 
parameters so that the element error e c - might carry 
more weight in determining the total error J. These 
weighting factors could be used to normalize the 
errors instead of using but it is easier to weight 
the importance of elements relative to one another if 
they are already normalized. 

The primary reason for the multilevel, nonlinear 
optimization structure used in all the formulations 
being compared herein is to keep the aerodynamic 
dimension as small as possible without adversely af- 
fecting the overall system analyses. The next three 
subsections present outlines of the optimization pro- 
cesses used by each of the extended approaches. The 
overall optimization methodology for each extended 
approximation formulation is depicted in figure 3. 
Differences in the actual approximation formulations 
result in variations in the optimization process of 
selected coefficients, but all three use similar tech- 
niques. Each requires a multilevel optimization pro- 
cess in the sense that each contains a closed-form lin- 
ear least-squares solution for the free linear param- 
eters, which is performed inside an iterative search 
for parameters that enter the problem in a non- 
linear fashion. Furthermore, each method requires 
that the characteristic roots of the system matrix 
corresponding to the aerodynamic states be stable. 
This constraint is enforced during optimization of the 
nonlinear parameters in the aerodynamic model. The 
choice of method (ELS, EMMP, or EMS) for a par- 
ticular application depends on aspects such as the 
goodness of fit required, computer size and time 
available, and computer costs. 
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The optimization of the nonlinear parameters has 
been included in obtaining “best” approximations 
using each of the extended LS, MMP, and MS 
formulations of section 3. The order of fit (n^ for 
ELS; til-, j = 1 , . . . , + n g , for EMMP; or N for 

EMS) is initially selected a priori. Best fits are then 
determined using this number. The order of fit is 
then increased or decreased based upon the suitabil- 
ity of the current best fits, until more suitable ones 
are obtained. 

The optimization of the linear coefficients can be 
considered a first level of optimization. Because of 
the nature of the objective function and the form 
of the approximating functions, the e c - can be opti- 
mized with respect to the linear coefficients in closed 
form, eliminating the need for iterative solutions typ- 
ical in nonlinear optimization. The optimization of 
the nonlinear parameters may be considered a second 
level of optimization. It is necessarily an iterative 
method which employs an algorithm to optimally se- 
lect the nonlinear parameters. Each selection of non- 
linear parameters requires one linear optimization 
for ELS and EMMP, and more for EMS, to be per- 
formed to determine the optimum linear coefficients 
corresponding to the particular selection of nonlinear 
coefficients. 

Figure 3(a) summarizes the overall optimization 
processes for the ELS approach, which is column- 
independent; figure 3(b) summarizes the EMMP ap- 
proach, which is column-dependent; and figure 3(c) 
summarizes the EMS process, which incorporates 
two linear iteration loops within the nonlinear opti- 
mization. Initially, in all three approaches, the form 
of the RFA to be used, the columns (or column) over 
which the optimization is to be performed, conver- 
gence criteria, etc., are identified. The process is then 
started by selecting an initial set of lag coefficients 
{bi } o (or {b^} o) and design weights The vec- 

tors of coefficients {A m } z y which. minimize the linear 
equality-constrained objective functions e c -- are then 
determined for each element in the desired column 
or columns of the matrix Q. The nonlinear objective 
function J for the ELS and EMS approaches or Jj for 
the EMMP approach is evaluated and a new set of 
values for the design variables is determined by the 
nonlinear optimizer. The best set of (A m ) z -y’s are 
then recomputed using matrix techniques. For the 
EMS approach, the overall process to determine the 
b% is the same as for the ELS approach except for the 
iteration loop to determine the D and E which de- 
fines a converged minimum value for J with respect 
to the current nonlinear parameters. 

4.3.1. Extended least-squares (column-independent) 
optimization process. This section outlines the steps 


in the overall optimization process for the column- 
independent ELS approach of reference 8. In step 1, 
the design variables for the linear and nonlinear 
optimization are separated by considering as design 
parameters the nonlinear b £ separately from the lin- 
ear coefficients (A m )ij in equation (3.1). Figure 3(a) 
is the flow diagram of the optimization procedure for 
this ELS approach. 

Step 1: For a specified set { b i : bf> > 0}, perform 
linear optimization . 

Step la: 

Determine all the (A m )ij and A y, for each 

which minimize each constrained least-squares 

error and satisfy specified constraints: 

Sn IQijd^n) QijU^n) P 

8i J “ ~m7. 

+ 2Xfj[cj{A m } i j-C° ij 

where 

M ij = max{l, \Qij(ik n )\ 2 } 

QijiP ) = (^o)ij (-^l )ijP 

n L 

+ {M )ijP 2 + Y^ A i +2 )ij-—^ 

Ste p lb: 

Calculate the total objective function 

3 = ( 4 - 8 ) 
V i 1 

Step 2: Select new set { b : b^ > 0} usina the 
nonlinear optimizer. 

Repeat steps 1 and 2 until an optimum solution 
for J is reached. 

4.3.2. Extended modified matrix-Pade (column- 
dependent) optimization process. This section outlines 
the steps in the overall optimization process for the 
column-dependent EMMP approach. Figure 3(b) de- 
picts the optimization procedure for this EMMP ap- 
proach. The distinguishing difference between this 
optimization process and the column-independent 
case of section 4.3.1 is that the b£ are no longer the 
same for all elements in the aerodynamic force ma- 
trix Q. In this column-dependent case, they depend 
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a given column and are necessarily the same only 
for the elements of that column, rather than for 
all the columns. This difference in formulation is 
accomplished simply by a redefinition of the non- 
linear objective function, which is now a column 
objective function identified as Jj. The optimization 
is now performed as a sequence of linear problems 
over a specified column for a fixed set of b ^ which 
are optimally selected by the nonlinear optimizer. 

Stev 1: For a specified set and a set {b^. : b^. > 
0}, verform linear optimization. 


Stev la: 


Determine all the (A m )ij and A z -y, for each ele- 
ment in column j of Q, which minimize each con- 
strained least-squares error and satisfy specified 
constraints: 


e 


c 

ij 


^2 n I Q ij ) Q jj [}k n ) | 2 

M{j 


+ 2A 




Cj 



where 


Mij = max{l, \Qij{ik n )\ 2 } 


QijiP ) — + (-^1 )ijP 



Stev lb: 

Calculate the column objective function 
J 3 = W ij £ ij 

Step 2: Select new set {6^. : b^ > 0} using the 
nonlinear oytimizer , 

Reveat steps 1 and 2 until an optimum solution and 
ni- for column j have been found. 

Reveat entire process for each column j — 1,..., 
+ n g . 


433 . Extended minimum-state optimization pro- 
cess . This section outlines the steps in the overall op- 
timization process for the EMS approach (fig. 3(c)). 


In the minimum-state formulation, the numerator co- 
efficients of the lag terms are the product elements 
of D and E, namely, 

(^£+2 )ij = D^Ee. (For t = 1, . . . , N) (4.9) 

Therefore, there are essentially two sets of linear 
coefficients in each minimum-state approximation, 
and the minimum-state optimization is necessarily 
a three-step process. First, an initial set {fy}o of el- 
ements for R and design weights {wij} are selected 
and the matrix D is set equal to an initializing ma- 
trix whose rank is the minimum of and -j- n g . 
Step la in the optimization process is a linear con- 
strained least-squares optimization which determines 
the best E, given D. Using this E in step lb, an- 
other linear constrained least-squares optimization is 
performed to determine Aq, Ai, A 2 , and a new D. 
Thus steps la and lb express the separation of the 
design variables for the two linear steps in the overall 
process for the minimum-state optimization. Step lc 
computes the total approximation error J for the re- 
sults of steps la and lb. Steps la to lc are repeated 
until the overall objective function J converges using 
the current set of b £ or hits an alternative stopping 
criterion, such as a maximum number of iterations. 
Step 2 then uses the nonlinear optimizer to select 
a new set of b £. The entire process continues until 
an optimal value for the overall objective function is 
reached. 


Step 1: For a specified set {bg : b £ > 0}, perform 

linear optimization. 

Step la: 

For j = 1 , ..., + n g (using the initial 

or previously determined D), determine all the 
(A m ) z y, Ely and A z y, which minimize the con- 
strained least-squares error e c - simultaneously, 
for all z = l, . . . , where 

Qij{p) = {Ao)ij + (Ai )ijP 

+ {M )»;p 2 E {j+bl) ( 4 - 10a ) 


Stej) lb: 

For i = 1, ri£ (using previously determined 
E), determine all the {A m )ij, A'£> and \ji which 
minimize the constrained least-squares error £? • 
simultaneously, for all j = 1, n^ + n g where 
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Qij(p) — (Ao)ij + {Al)ijP + ( A2)ijP 2 

^ ( E lP \ 

+ £ D " (?+*;) (4,10b> 

Stev lc: 

Compute the total objective function 
J = 

V 3 * 

Repeat Steps la-lc until J converges, using the cur- 
rent {b e }. 

Steg 2j Select new set {be : be > 0} using the 
nordinear ogUmizer which reduces the total objective 

j 'unctwJL £ ■ 

Repeat steps 1 and 2 until an optimum solution for 
J is reached. 

Figure 3(c) summarizes the overall process which in- 
cludes the dual-level linear optimizations for E, D, 
and the A m . The additional shaded area in fig- 
ure 3(c) indicates the logic required by the minimum- 
state optimization process. 

4.4. Side Constraints for Nonlinear 
Optimization 

Since the be (or be-) are allowed to vary, it is 
necessary to impose side constraints on them. These 
lag coefficients must be greater than 0 in order to 
ensure system stability as a result of introducing the 
related aerodynamic states into the state equations. 
Also, it is frequently desired to restrict the range of 
variation to that in the neighborhood of the range 
of frequencies over which tabular data are available; 
that is, 

0 < Le <be <U e (For l = 1, . . . , n-i) (4.11a) 

or, for column-dependent optimization 

o < L( j < be j < Ue j (For j = 1, . . . , n^ + n g , 

i=\, ..., n L .) 

(4.11b) 

These side constraints are enforced in one of two 
ways. The first is by way of an inverse sinusoidal 
transformation of the design space [Le, Ue\ defined by 
equation (4.11a) onto the real line segment [— 1,+1] 
(ref. 26). The relationship between the two line 
segments is 


+ ^ (4.12a) 

where 

-1 < Zfi < 1 (4.12b) 

As can be observed, the constraint boundaries are 
mapped onto z = ±1. This transformation ensures 
that the side constraints are always satisfied. How- 
ever, the restrictions on the defined by the in- 
equality (4.12b), are not strictly enforced by the non- 
linear optimizer. Consequently, an occasional oscilla- 
tion problem between successive values of the design 
parameters arises because of the multivalued charac- 
teristic of the transformation. Various methods could 
be applied to avoid this problem; since it rarely oc- 
curs, however, a penalty function formulation of the 
nonlinear objective function (refs. 13, 14, 27, and 28) 
is employed to enforce the side constraints when it 
does. 

Two examples of available penalty functions are a 
wall function and the extended-interior penalty func- 
tion. By wall, we mean that the function essentially 
hits a wall when constraints are violated. This tech- 
nique is implemented by defining the objective func- 
tion as very large relative to its normal range of 
values. This has extreme discontinuities which can 
sometimes cause convergence problems and usually 
requires that an optimization process start within the 
feasible design space. 

The second type of penalty function, proposed by 
Haftka and Starnes (ref. 27), is an extended-interior 
formulation. It places no initial restrictions and does 
not inject discontinuities into the convergence pro- 
cess. This is normally the preferred type of penalty 
function if convergence problems do arise from the 
constraint violations since convergence in this case is 
smoother. The sinusoidal transformation and both 
the wall and the extended-interior penalty function 
approaches to constrain the nonlinear parameters 
have been used in the optimizations presented here. 

5. Numerical Application to the DAST 
ARW-2 

The three formulations described in this paper 
were all applied to a mathematical model of the 
DAST (Drones for Aerodynamic and Structural Test- 
ing) ARW-2 (Aeroelastic Research Wing Number 2) 
aircraft. The ARW-2 vehicle was designed to have 
a maximum gross weight of approximately 2500 lb 
and a wing with a supercritical airfoil having a span 
of 18.98 ft, an aspect ratio of 10.3, and a sweep an- 
gle of 25° at the quarter-chord. The wing was pur- 
posely designed to require flutter suppression, ma- 
neuver load alleviation, gust load alleviation, and 
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static stability augmentation in some regions of its 
flight envelope. 

5.1. Description of Mathematical Model 

A modal characterization of the ARW-2 was em- 
ployed which was obtained from a free-free vibra- 
tion analysis of the ARW-2. Twelve modes were re- 
tained for the symmetric degrees of freedom. These 
are plunge, pitch, and 10 symmetric elastic modes. 
Table 1 lists the natural frequencies associated with 
each of the in-vacuum elastic modes. For the results 
presented here, seven modes were employed. They 
were composed of 5 of the original 12 modes (-plunge, 
pitch, and elastic modes 1, 4, and 6), a rigid control 
surface rotation mode (stabilizer), and a sinusoidal 
gust mode. The deleted modes had little effect on 
flutter or rigid body stability and response charac- 
teristics. The actuator driving the control mode was 
considered reversible. 

A doublet lattice method (refs. 17 and 18) that is 
contained in the ISAC system (ref. 15) was used to 
obtain tabular values of the generalized aerodynamic 
force coefficients at a Mach number of 0.86 for the 
following set of reduced frequencies: 


Baseline- 1 (Bl): No equality constraints imposed, 
no nonlinear optimization 

a. 1 lag coefficient (6 states) 

b. 2 lag coefficients (12 states) 

c. 3 lag coefficients (18 states) 

d. 4 lag coefficients (24 states) 

Figure 5(a) shows the approximations with 4, 3, 
2, and 1 lag term for a single (typical) aerodynamic 
element (lift due to the fourth elastic mode) with lag 
coefficients evenly spread over the reduced-frequency 
range. The open circles are the actual tabular values, 
and the solid circles are the corresponding values of 
the approximation at the same values of frequency. 
The fits deteriorate with a reduction in the order of 
fit. Other aerodynamic elements (not shown) indi- 
cate a similar trend. For each of these unconstrained 
cases, the total approximation error J for all aero- 
dynamic elements is shown in figure 6. 

5.2.2. Total approximation errors for constrained , 
unoptimized approximations . The second baseline 
case again uses the ELS formulation for different 
numbers of lag coefficients in the approximation, but 
with the following equality constraints imposed: 


{k n } = {0,0.005,0.01,0.05,0.1,0.2,0.3,0.4,0.5, 

0.6, 0.8, 1.0} 


1. The steady-state values (at k = 0) are con- 
strained to agree for all the coefficients: 


and uj n = ( 2ufc)k n where c = 1.956 ft and 
u = 874.75 ft/sec for the cases studied. The panel- 
ing of the lifting surfaces is shown in figure 4(a). The 
circles indicate points at which mode shape data are 
defined. 

Frequency responses of pitch due to control sur- 
face rotation and c.g. acceleration due to control sur- 
face rotation were calculated using ISAC and are pre- 
sented later in this section to illustrate the accuracy 
of the modeling. Figure 4(b) depicts the locations of 
the control surface (stabilizer) and the sensors (pitch- 
rate gyro and c.g. accelerometer) used in these fre- 
quency response analyses. 

5.2. Baseline Comparisons 

As a way of introducing the effects on the RFA’s 
of (1) the number of lag terms and (2) the presence 
or absence of constraints, some baseline comparisons 
are presented using the ELS formulation. 


Qij{ 0 ) = Qij (0) (For i = 1 , . . . , 

and j = 1, . . . , n^ + rig) 

(5.1) 

2. The steady-state slopes (at k = 0) are con- 
strained for the force coefficients due to -plunge 
(mode 1) and pitch (mode 2): 


dQilip) _ Qj 2 (0) 

dp p =o b9 n 


(i 9 n = -0.00624865) 
(5.2) 


dQilip) 

dp 


p=0 


_ Im[Q t2 (ffr 2 )] 

k 2 


(5.3) 


where assumptions are that /ci = 0 and /c 2 is small. 

3. For the aerodynamic force coefficients due to 
the flexible modes, control mode and gust velocity 
are all constrained to agree with interpolated tabular 
data at the reduced frequency ( kf = 0.127) near 
flutter: 


5.2 .1. Total approximation errors for unconstrained , 
unoptimized approximations . The first baseline case 
uses the ELS form for approximating the unsteady 
aerodynamic forces; different numbers of lag coeffi- 
cients are used in the approximation, and no con- 
straints are imposed: 


QijUkf) — Qijiikf) (For i — 1, . . . , ri£ 

and j = 3, . . . , ri£ + n g ) 
(5.4) 

The rationale for constraint 2 is discussed in refer- 
ence 10 and for convenience to the reader is included 
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in appendix B. Thus the second baseline case is as 
follows: 

Baseline-2 (B2): Equality constraints imposed, 
no nonlinear optimization 

a. 1 lag coefficient, constrained (6 states) 

b. 2 lag coefficients, constrained (12 states) 

c. 3 lag coefficients, constrained (18 states) 

d. 4 lag coefficients, constrained (24 states) 

Figure 5(b) shows a constrained approximation 
for the same element of Q as shown in figure 5(a), 
but with constraints described as for equation (B2). 
In each case the fit at and near k = 0 is very good 
and each of the approximating curves also crosses the 
interpolated curve at kf. 

5.2.5. Comparison of constrained and uncon- 
strained results . Comparing figure 5(a) with 5(b) 
shows that there is a significant difference between 
the unconstrained and constrained fits for higher fre- 
quencies, and the fits deteriorate drastically for the 
constrained cases as the order of fit decreases. Both 
parts of figure 5 demonstrate that the approximating 
curve can be improved by increasing the number of 
b^. Although not shown, the fit can also be improved 
by reducing the frequency range over which the fits 
are required, but this reduces the robustness of the 
approximation. Reduction in frequency range for a 
column may be desirable if physical considerations 
dictate that the mode has negligible effect outside 
some frequency band. Figure 5 also indicates that 


5.5.1. Total approximation errors for constrained , 
optimized approximations . Table 3 and figure 7 

compare the three extended approximation methods 
solely in terms of J defined by equations (4.7), for 
each of the optimized, constrained cases as a func- 
tion of the aerodynamic dimension. Typically, when 
optimizing this error function having weighting fac- 
tors wij , one would weight important elements more 
heavily relative to less important ones and perform 


the number of lag coefficients used in the approxima- 
tions can be a critical factor in the analyses for which 
these approximations are being determined. 

Table 2 lists the total approximation error J 
with respect to the aerodynamic dimension for both 
baseline cases. The approximation error increases 
with the imposed constraints by almost a factor of 11 
for the 6-state case, and a factor of 1.4 for the 24-state 
case. Figure 6 is a plot of the data in table 2 
and indicates more clearly that optimization of the 
nonlinear coefficients is indeed desirable, especially 
when constraints are imposed. It is less necessary as 
the aerodynamic dimension is increased. 

5.3. Optimized Comparisons 

The optimizations of the nonlinear coefficients 
using the three rational function formulations de- 
scribed in this paper were performed with the equal- 
ity constraints imposed as in the baseline-2 cases, 
resulting in the following optimized cases to be 
compared: 

ELS: Equality constraints (eqs. (5.1)— (5.4)) 

imposed 

EMMP: Equality constraints imposed, the same 

order of fit assumed for each column 

EMMP/v: Equality constraints imposed, the num- 
ber of denominator coefficients varied 
for each column 

EMS: Equality constraints imposed 


the fits only over the frequency range pertinent to 
the problem being analyzed. The weighting factors 
in the error function for this example, however, were 
chosen to be 1.0, and the fits were done over the 
entire frequency range for all elements. 

Note that there are two columns in table 3 
and two curves in figure 7 that correspond to the 
EMMP formulation. For the column and curve 
labeled “EMMP,” each column has the same 


Cases Compared 


ELS 

EMMP 

EMMP/v 1 

EMS 

6 states (1 lag) 

12 states (2 lags) 
18 states (3 lags) 
24 states (4 lags) 

7 states (1 lag) 
14 states (2 lags) 
21 states (3 lags) 
28 states (4 lags) 

7-28 states 

(varying number 
of lags per 
column) 

4 states (4 lags) 

6 states (6 lags) 

7 states (7 lags) 

8 states (8 lags) 

9 states (9 lags) 

10 states (10 lags) 
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number n^. of lag coefficients 6^., whereas for the 
column and curve labeled “EMMP/v,” n L . differs 

between columns. The EMMP/v case, which em- 
ploys the full flexibility of the matrix-Pade approx- 
imation approach, is superior to the EMMP case 
where n^. is the same for all columns. As expected, 
the ELS method, which is required to have the same 
b£ for each element, has larger approximation er- 
ror, although slight in this case, than either of the 
EMMP approaches for a given aerodynamic dimen- 
sion, and the EMS approach achieves the lowest er- 
ror of all for a given aerodynamic dimension. One 
can see by comparing figure 6 with figure 7 or by 
looking at table 3 that for each formulation, the to- 
tal error is improved over the corresponding base- 
line case. Figure 7 also shows that for each for- 
mulation, the errors decrease as the aerodynamic 
dimension increases, but for sufficiently large aero- 
dynamic dimension, the difference in error between 
the approaches becomes insignificant. This is espe- 
cially true of the ELS and EMMP approaches. 

Another result, which corroborates findings of 
reference 7, can be seen by comparing, for each 
formulation, the aerodynamic dimension required to 
achieve a given level of error. It is shown in figure 7 
and table 3 that the fit errors for the minimum- 
state formulation for low aerodynamic dimensions 
are comparable to the fit errors for the EMMP and 
ELS formulations at significantly higher aerodynamic 
dimensions. For example, the 7-state EMS error 
falls right between the 12-state ELS error and the 
14-state EMMP error, and the 10-state EMS error 
is approximately equal to the EMMP/v-17 error and 
only about 12 percent higher than the 24-state ELS 
error. 

The ranking from left to right of the four curves 
in order of total error for a given aerodynamic 
dimension seen in figure 7 will hold in general, but 
the relative spacing between curves will be problem 
dependent. In general, the ELS would tend to lie 
further to the right of the EMMP curves. The total 
approximation error, however, is only one measure of 
the goodness of the RFA’s. The approximations for 
individual elements must also be examined. Three of 
the optimized cases, namely, 

ELS with 4 lag coefficients, 24 states (ELS-24) 

EMMP with 3 lag coefficients per column, 

21 states (EMMP-21) 

EMS with 10 lag coefficients and 10 states 
(EMS-10) 

were selected for further comparisons since the 
EMMP-21 and the EMS-10 have nearly the same 


total approximation errors, and the ELS-24 has an 
error which is closest in value to the other two. 

Cases where one aerodynamic element is difficult 
to fit but too important to weigh less than the others 
can occur. This can complicate the process of obtain- 
ing satisfactory RFA’s by disproportionally influenc- 
ing the selection of the nonlinear parameters, and for 
the minimum-state, even the linear parameters. This 
is illustrated in the next section. 

5.3.2. Comparison of fits for selected GAF ele- 
ments. Figure 8 shows the lift and pitching moment 
due to each of the seven modes selected for these 
studies (-plunge, pitch, three elastic modes, control 
surface rotation, and sinusoidal gust). In figure 8, 
the short-dashed lines are interpolated fits through 
the tabular data, indicated by circles, and represent 
the “truth model.” The other lines are the approxi- 
mating rational functions, and the solid symbols cor- 
respond to the values of the approximating functions 
at the same frequencies as the tabular data. As can 
be seen from these samples, some of the elements 
are much harder to fit than others. Since the selec- 
tion of the &£ is based on the approximation errors 
for all elements, the ELS and EMS approaches are 
more sensitive to the hard-to-fit elements than the 
EMMP approach. The hard-to-fit elements affect the 
selection of b £ only for a particular column in the 
EMMP approach. The EMS is the most sensitive 
since the numerator coefficients for the lag terms are 
coupled to all the columns through the matrices D 
and E. The driving element, shown in figure 8(n), for 
the approximations being compared is the pitching 
moment due to gust. For the set of constraints se- 
lected, the EMS actually gives the best fit for this 
element in the low frequency range, but in doing so, 
the fits for the other elements are severely degraded. 

Table 4 also demonstrates this phenomenon. By 
comparing the first three columns showing the col- 
umn error Jj for each of these three cases (ELS-24, 
EMMP-21, and EMS-10), it can be seen that the 
EMS-10 case actually obtains the least error for the 
gust column (column 7), but it does so at the expense 
of all the other columns. Results of the EMMP/v-17 
case (only one lag coefficient for column 2, two lag 
coefficients for columns 3, 5, and 6, three lag coeffi- 
cients for columns 1 and 4, and four lag coefficients 
for the gust column) are shown in the fourth col- 
umn of table 4 and in figure 9. Comparing these 
with the EMMP-21 case with three lags per column 
shows slightly worse fits for the forces due to pitch 
and elastic modes 1 and 6 and the control mode, 
which have fewer lag coefficients, the same for those 
due to plunge and elastic mode 4, which have the 
same lag coefficients, but are better for those due to 
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gust, which has one more lag coefficient. The over- 
all error is improved by 7 percent with a reduction 
in aerodynamic dimension of 19 percent. This can 
be explained by the fact that fewer lag coefficients 
are being used for the easier-to-fit elements in order 
to reduce the aerodynamic dimension without signif- 
icantly affecting the total approximation error. The 
additional lag used for the forces due to gust helps to 
significantly reduce the error in fit for these with only 
one additional state. Again, however, the gust col- 
umn even in these cases (EMMP-21 and EMMP/v- 
17) is driving the total approximation error. Note 
also the fits overall for the EMMP/v-17 are about 
as good as the EMS-10 fits as indicated by the fact 
that the total errors (2.409 and 2.420, respectively, 
in table 3) are nearly the same. 

By relaxing the constraints on the gust column 
so that the elements are not constrained to fit at the 
flutter frequency, the curves for the other elements 
can be improved considerably. Since the EMS-10 was 
the most sensitive to this constraint, optimizations 
were performed for this approach with the flutter 
constraint removed for the gust column as well as the 
control column. The results are shown in figure 10. 
In general, the approximations shown are markedly 
improved. Removal of the constraint to match gener- 
alized aerodynamic forces at the flutter frequency for 
these two columns is reasonable since the stabilizer 
will be applied only to control low-frequency rigid 
body motion; and for points removed from flutter, 
loads due to gust have most of their power concen- 
trated within the rigid body frequency region. In 
fact, one could also reduce the frequency range of 
interest, provided that at least as many equations 
are retained as there are independent variables. For 
the higher order EMS cases, however, this reduction 
in frequency range adversely affects the least-squares 
linear solution. 

5.3.3 . Comparison of selected frequency responses. 

Other measures of goodness of fit, such as frequency 
responses, may also be used to demonstrate the ad- 
equacy of candidate p-plane approximations. To 
illustrate the effects of the three RFA’s on “final an- 
swers” (i.e., solutions of the aeroelastic equations of 
motion), each of the three RFA’s were implemented 
in the aeroelastic equations of motion. Pitch and c.g. 
acceleration frequency response functions due to con- 
trol deflection were computed for the following flight 
condition: 

Mach number = 0.86 

Velocity, u = 874.75 ft/sec 

Dynamic pressure, q = 411.5 lb/ft 2 

Figures 11 and 12 contain these frequency response 


functions. Figure 11 shows the pitch frequency 
response function and figure 12 shows the c.g. acceler- 
ation response due to control deflection for a selected 
flight condition. Shown on each figure are a baseline 
case for which no aerodynamic force approximations 
using rational functions are employed, ELS-24, which 
had the least total error of the three p-plane approx- 
imation cases, and EMMP-21 and EMS-10, both of 
which had nearly the same value for the total approx- 
imation error. There is good agreement between the 
frequency responses for all three cases, but the ELS 
is not the best, contrary to what might be expected. 

There may be several solutions for the b £ which 
will give similar total approximation errors. As 
implied by the results in figures 11 and 12, it may be 
desirable to consider other measures of how good the 
approximations are than just the total approximation 
error itself. The accuracy of key frequency response 
functions for various fits is another way to assess the 
adequacy of candidate p-plane approximations. 

In summary, these results illustrate some of the 
factors that should be considered and indicate that 
considerable engineering judgment may be required 
to obtain both sufficiently accurate RFA’s and low 
aerodynamic dimension. 

5.4. Hidden Costs and Related Problems 

The advantage of optimizing the nonlinear coef- 
ficients in each of the formulations presented in or- 
der to reduce errors for a specified aerodynamic di- 
mension is clear. The ability to significantly reduce 
the aerodynamic dimension itself using the formula- 
tion of reference 8 is especially promising. However, 
these advantages do not come without cost. A hidden 
cost here, of course, is the increased computer cost 
needed to perform the optimizations. The computer 
costs in the cases compared here are approximately 
as follows: 

ELS = 20 times baseline cost 

EMMP = 80 to 100 times baseline cost 

EMMP/v = 700 to 1000 times baseline cost 

EMS = 2500 to 3000 times baseline cost 

These costs may ultimately be offset because the 
resulting matrices used later in aeroelastic and aero- 
servoelastic analyses are smaller. 

Another hidden cost, in the minimum-state opti- 
mization process, is an increased likelihood of compu- 
tational limitations associated with the EMS due to 
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the size of the constrained least-squares linear system 
of equations to solve. These limitations can cause 
convergence problems in the dual linear optimization 
used to calculate the total error, which lead to prob- 
lems within the nonlinear optimization since the ap- 
proximation error J has not converged. This did not 
occur, however, for the cases presented. 

The desirability of enforcing constraints is based 
on the physical aspects of the problem being consid- 
ered. In most instances, flexibility in selecting con- 
straints results in more realistic models and more 
reliable analyses. Since the larger costs for the 
minimum-state formulations are due primarily to 
the internal iterations in the dual linear optimiza- 
tions, fixing the constraints for all elements as Karpel 
(ref. 8) did would reduce some of these costs. He en- 
forced the same three constraints for all columns, and 
using these constraint equations, he solved for the 
matrices Aq, Aj, and A 2 explicitly, thereby reduc- 
ing the number of unknowns in the linear optimiza- 
tion. Use of a similar technique on a per-column basis 
to include some of the selected equality constraints 
provides a mechanism to reduce the number of un- 
knowns in the linear optimization while retaining the 
advantages of being able to select desired constraints. 

6. Concluding Remarks 

Three currently used formulations (namely, a con- 
ventional “least-squares” (LS), a modified matrix- 
Pade (MMP), and a “minimum-state” (MS)) for 
approximating unsteady aerodynamic forces in the 
equations of motion of a flexible aircraft with rational 
functions in the Laplace domain are reviewed. Op- 
timization schemes have been applied to determine 
the parameters in each formulation. Furthermore, 
a number of selectable constraints on the approxi- 
mating functions have been included. This paper ex- 
tends the modified matrix-Pade and minimum-state 
approximation methods to the same level of con- 
straint selection and uses the same nonlinear op- 
timization techniques as the extended least-squares 
method previously developed by the authors. The 
three extended methods (ELS, EMMP, and EMS) 
were applied to an aeroservoelastic research model to 
provide comparative evaluations. Results from the 
applications show that the number of aerodynamic 
states employed in representing the equations of 


motion in linear time-invariant form can be signif- 
icantly reduced from the number required by the 
conventional constrained least-squares approach. In 
fact, the EMMP approach and the EMS approach 
resulted, respectively, in 24 percent and 67 percent 
fewer aerodynamic states than a corresponding base- 
line constrained least-squares approach having a sim- 
ilar total approximation error in which nonlinear pa- 
rameters were selected a priori. Increases in com- 
puter costs as a result of the optimizations performed 
for ELS, EMMP, and EMS were approximately 20, 
700, and 2500 times the cost of the unoptimized base- 
line approach. In summary, it has been shown that 

1. Optimization of the nonlinear parameters 
improves the approximations in all three 
formulations. 

2. Use of the full flexibility in the modified matrix- 
Pade formulation (i.e., variable number of lag 
terms per column) significantly reduces the aero- 
dynamic dimension of the state-space equations 
of motion from that required when each column 
has the same number of lag terms. 

3. For similar total approximation errors, the aero- 
dynamic dimension for EMMP is less than that 
for ELS. However, if control surface actuators can 
be approximated as irreversible rather than re- 
versible, the aerodynamic dimension for ELS de- 
creases and the difference between the two is less. 

4. The minimum-state formulation shows great 
promise in allowing a significant reduction in the 
aerodynamic dimension of the state-space equa- 
tions of motion. The approximation errors for 
the minimum-state formulation for low aerody- 
namic dimensions are comparable to those for the 
other two formulations with higher aerodynamic 
dimension. 

Regardless of the formulation used, to obtain suf- 
ficiently accurate rational function approximations 
with low aerodynamic dimension, optimization has to 
be coupled with engineering judgment to determine 
the relative importance of individual aerodynamic 
elements, critical frequency ranges, and appropriate 
constraints to impose. 

NASA Langley Research Center 
Hampton, Virginia 23665-5225 
March 10, 1988 
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Appendix A 

Nonlinear, Nongradient Optimizer 

A.l. Sequential Simplex 

The algorithm used to perform the nonlinear op- 
timization in all the RFA formulations presented 
in this paper is a sequential simplex (or polytope) 
method developed by Nelder and Mead (refs. 12-14). 
It has been used in control system design for flex- 
ible aircraft with active controls (refs. 20 and 21), 
along with gradient methods such as the CONMIN 
feasible direction method (ref. 22), the Davidon- 
Fletcher-Powell variable metric (refs. 23 and 24), and 
the newer Davidon optimally conditioned method 
(ref. 29). 

The nongradient Nelder-Mead sequential simplex 
algorithm is reliable and efficient. The method is 
simple to use and robust in its ability to handle the 
nonlinear optimization problem to which it has been 
applied. Its adaptive nature in moving away from 
“high” points requires minimal effort on the part of 
the user to “fine tune” the problem. It is a more sta- 
ble and robust algorithm, numerically, than some of 
the gradient-based algorithms to which it has been 
compared. Since computation of gradients using fi- 
nite differencing can be costly in large applications 
programs for which closed-form gradients are either 
not available or not practical, this nongradient al- 
gorithm has proven to be invaluable. It lacks the 
initial sensitivities of gradient-based methods. But, 
the final two low values and the step size between 


them do provide some measure, although somewhat 
obscure, of sensitivity at convergence. At this point, 
the gradient could be computed if desired. 

The basic mathematical justification of conver- 
gence for this algorithm is based only on the conver- 
gence of a monotonically decreasing series. Although 
its order of convergence depends on the function be- 
ing optimized, it has been found to converge rapidly 
in numerical applications. 

A.2. Description of Algorithm 

Figure 13 depicts the algorithm for a two- 
dimensional design space. The algorithm starts with 
a simplex of points in the design space (A ABC). The 
objective function is evaluated at each vertex in the 
simplex and the highest valued point is identified (A). 
A line of projection through the centroid of the oppo- 
site side (namely, the other n vertices) is determined 
and the objective function is then evaluated at the 
reflected point (E) on the projection line. Depend- 
ing on the relative value of the function at this point 
and the other points in the simplex, a single exten- 
sion (F), a retention (E), an exterior contraction (H), 
or an interior contraction point (G) is then identified 
at which the objective function is evaluated and a 
new simplex is determined. The actual decision pro- 
cess is detailed in references 12-14, although the code 
as listed (ref. 13) deviates from the decision process 
as described in the reference with respect to those 
steps taken when equalities hold. 
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Appendix B 

Rigid Body Slope Constraint 

Linear relationships exist between certain func- 
tions of some generalized force elements at steady 
state (k = 0) which can be included as constraints 
on the rational function approximations; see refer- 
ence 10. An example is illustrated here for the 
-plunge and pitch modes. As indicated previously 
(eq. (2.4)) the generalized forces in the reduced- 
frequency domain are 

^ . f f A PAx, y, s) 

Qij(s) = jj — • — - Zi{x, y) dS (Bl) 


which can also be expressed as 





J k 1 {x,y y x , ,y',k)- ( X ^ dS' 


dS 


(B2) 


where k is an influence function that defines the 
downwash at [x* ,y f ) due to a unit pressure distur- 
bance at {x,y). The deflections (z z , zq) and normal- 
ized downwash (w z ,wq) for the -plunge and pitch 
modes oscillating at frequency w, can be expressed 
in terms of the mode shapes (z\, z^) and their slopes, 
respectively, as follows: 


—plunge (mode 1) 

z z = ziti = = e ( iu / b ) kt ' 

vo z _ 10J ( t - W £) _ ik 

— e — ~r z z 
u u b 


Taking the Fourier transform, 

z z (k') = 6(k' - k ) 

zy(k ) = ( x Xcg)9n,6 (k fc) 

— (k 1 ) = -r-S(k' - k) 
u b 

X ^{k') = jz e {k')- 0 n 6 {k'-k) 


(B4) 


Substituting these into equation (B2) yields the un- 
steady aerodynamic forces in the reduced-frequency 
domain: 


Qn{k) = J Zi(x,y) [/ K 1 (x,y,x' ,y' ,k) y dS 1 
= J z i j J K~ 1 (x,y,x',y',k)(ik) dS 1 


dS 

dS 


Hence, 

dQn(k) 


dik 

Also, 


c =o = G) / Zi \j k 1 ( x <y^ x '>y'< 0 ) dS ' dS 


QMk) = J Zi(x,y ) 

x [/ K~ 1 {x,y,x\y\k) - 9n^J dS' dS 

Hence, 

Qt 2 ( 0 ) = -On J Zi(x,y) J K _1 (x,7/,x , ,j/ , ,0) dS 1 dS 


pitch (mode 2) 


> (B3) 


Z0 - zifr - (x - x C g)O n e lut 
= -(x-x cg )0„e (i “ /,,)fct 

u J-=™z e - 9 n e iut 
u u 

= * zg - e n e {iu/b ^ kt 
0 


Therefore, 


9Qji{k) 

dik 


k=Q 


= K qmo) 


A similar relationship exists between lateral displace- 
ment and sideslip or yaw. Explicit enforcement of 
these constraints improves the prediction of the low- 
frequency characteristic roots. 
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Symbols 

(Ao )»jj coefficients used in RFA for an 

(Ai)y , . . . element of matrix Q (see eq. (2.6)) 

Ao, Ai, ... matrix of coefficients used in RFA’s for 
matrix Q (see eq. (3.1)) 

{A m }ij column vector of coefficients 

Aq,Ax, ... for the (i, j) th element of 
Q (see eq. (4.3)) 

bfrbij constant coefficient in £th denominator 

of partial fraction expansion of RFA 
for element Q z y (p) , of matrix Q(p); 
the partial fraction term in the RFA 
containing this coefficient is commonly 
referred to as an “aerodynamic lag” 
term and the bp as “aerodynamic lag” 
coefficients or simply, “lag” coefficients 
(see eqs. (2.6) and (3.11)) 

c mean aerodynamic chord; reference 

length is c/2 

c.g. center of gravity 

Cj matrix multiplier of coefficients used 

in defining linear equality constraints 
for jth column of Q (see eq. (4.3)) 

C?- matrix of constants used in defining 

linear equality constraints for column 
of Q (see eq. (4.3)) 

D pre-multiplier numerator matrix 

for minimum-state formulation (see 
eqs. (3.17) and (3.24)) 

E post-multiplier numerator matrix 

for minimum-state formulation (see 
eqs. (3.17) and (3.24)) 

F matrix of generalized aerodynamic 

forces due to aircraft and control 
surface motions and due to gusts (see 
eq. (2.2)) 

Fhm x n 8 matrix of modal coefficients 

converting hinge-moment outputs to 
generalized forces (see eq. (2.3)) 

G matrix of damping coefficients used 

in Lagrangian formulation of the 
equations of motion (see eq. (2.2)) 

hj{p) denominator polynomial for jth 

modified matrix-Pade rational function 
approximation (see eqs. (3.12)) 

hj m mth coefficient of hj(p) (see eqs. (3.13)) 


H j state-space matrix multiplier of states 

in phase-canonical state-space real- 
ization of hj(p ) which defines the cor- 
responding augmenting aerodynamic 
state vector X a ^. (see eqs. (3.13)) 

I identity matrix whose dimension 

depends on equation in which it 
resides 

i complex variable, \/— T 

J total approximation error in approxi- 

mations for all elements of Q; specif- 
ically, weighted sum of square errors 
in approximations to elements of Q as 
a function of lag coefficients currently 
being used for approximations (see 
eq. (4.7b)) 

Jj total approximation error in approx- 

imations for all elements in jth col- 
umn of Q; specifically, weighted sum 
of square errors in approximations to 
elements in jth column of Q as a func- 
tion of lag coefficients currently being 
used for that column (see eq. (4.7a)) 

k nondimensionalized reduced frequency, 

u>c/2u 

k f reduced frequency corresponding to 

flutter 

k n reduced frequency at which general- 

ized forces are computed 

K generalized stiffness matrix used in 

the Lagrangian formulation of the 
equations of motion (see eqs. (2.2) 
and (2.3)) 

M generalized mass matrix used in 

the Lagrangian formulation of the 
equations of motion (see eqs. (2.2) 
and (2.3)) 

factor used to normalize element 
error for nonlinear optimization (see 
eq. (4.1b)) 

N number of lag terms and the aerody- 

namic dimension in the minimum- 
state formulation 

n a aerodynamic dimension 

n cj number of constraint equations for 

RFA of elements in ^th column of Q 
(see eq. (4.5b)) 

rig number of gust modes retained in 

equations of motion 
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nL 


n Lj 


n 6 

n 6 , 

n Z 


n £t 


P 

Q 

Q 

Q 

R 


s 

S 

t 

t hm 


V; 


Wg,Vg 


number of lag terms (equivalent to 
number of partial fractions and the 
order of the overall denominator 
polynomial) in least-squares RFA (see 
eqs. (2.6) and (3.1)) 

number of lag terms in matrix-Pade 
RFA’s for jth column of aerodynamic 
forces (see eq. (3.11)) 

number of actuators 

number of irreversible actuators 

number of modes retained in the 
equations of motion to define the 
motion of the vehicle, including 
control modes 

number of modes retained in the 
equations of motion to define the 
motion of the vehicle, excluding 
the control modes corresponding to 
irreversible actuators 

nondimensionalized Laplace variable, 
(c/2u)s 

dynamic pressure 

matrix of generalized force coefficients 
in Laplace domain (see eq. (2.4)) 

matrix of rational function approxima- 
tions to elements of Q (see eq. (2.7)) 

diagonal matrix of denominator 
coefficients used in the minimum- 
state formulation of the RFA’s (see 
eqs. (3.17)-(3.24)) 

complex Laplace variable 

total lifting surface of the aircraft (see 
eq. (2.4)) 

time variable 

vector of hinge moments output by the 
actuator (see eq. (2.3)) 

reference velocity 

state-space matrix multiplier of 
inputs in phase-canonical state- 
space realization of hj(p) which 
defines the corresponding augmenting 
aerodynamic state vector X a ^- (see 
eqs. (3.13)) 

vertical and lateral gust velocities 


Wij weight applied to the approximation 

error for the (i,j) th element of Q 
in determining the total error J or 
column error Jj (see eqs. (4.7)) 

W j matrix of polynomial coefficients as 

defined by equations (3.12), along with 
equation (3.14) used in matrix-Pade 
RFA formulation corresponding to jth 
column of Q 

x x-coordinate in spacial coordinate 

system 

X a aerodynamic state vector in aug- 

mented LTI state-space equations (see 
eqs. (3.4), (3.13), and (3.20)) 

y y-coordinate in spacial coordinate 

system 

z, z(x,y,i) ^-coordinate in spacial coordinate 

system representing the deflection at 
a point (x, y) at a given time t (see 
eq. (2.1)) 

Zi, Zi(x,y) time independent mode shape compo- 
nent of deflection 2 (see eq. (2.1)) 

ctg nondimensional gust velocities, 

[Wg/u Vg/u] T 

APj(x,y,s ) change in the lifting pressure due 
to changes in the jth generalized 
coordinate at a point (x, y ) on the 
lifting surface (see eq. (2.4)) 

£ij(ikn) approximation error between the 

tabular value and the approximation 
value of a generalized force at p = ik n 
(see fig. 1) 

unconstrained least-square error in 
fitting an element Qij of Q as given 
by equation (4.1a) 

constrained approximation error (in 
the Lagrangian formulation of the 
error) for of Q as defined by 
equation (4.4) 

vector of Laplace transformed time 
derivatives of generalized coordi- 
nates and gust inputs, as defined in 
equation (3.4) 

vector of inputs in state-space realiza- 
tion of least-squares RFA formulation, 
as defined in equation (3.4) 


V 

4 
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o 


Kj vector of size n C] of Lagrangian mul- 

tipliers in the Lagrangian formula- 
tion of the constrained objective func- 
tion for the (*,j)th element of Q (see 
eqs. (4.4)) 

w frequency of oscillation 

0 n normalizing angle in definition of rigid 

body pitch mode shape (see eq. (5.2) 
and appendix B) 

£i(t) generalized coordinate representing 

the amplitude of the corresponding ith 
component mode shape in the total 
deflection as a function of time (see 
eq. (2.1)) 

£* £, ( vector of generalized coordinates and 
its first two time derivatives used in 
Lagrangian formulation of equations of 
motion (see eq. (2.2)) 

£(s) vector of generalized coordinates in 

Laplace domain (see eq. (2.3)) 

vector of generalized coordinates ex- 
cluding control deflections correspond- 
ing to irreversible actuators 

Subscripts: 

g related to gusts 

i ith row of matrix Q 

j j th column of matrix Q 

t ^th partial fraction in RFA 


initial value 
6 related to controls 

f related to generalized coordinates 

Notation: 

[■] T transpose of a matrix [■] 

{■} column vector 

— bar over a function or matrix indicates 

that it corresponds to the dimensional 
Laplace s, rather than the non- 
dimensional p 

Acronyms: 

ELS extended least-squares 

EMMP extended modified matrix-Pade 

EMMP/v extended modified matrix-Pade with 
varying number of denominator 
coefficients per column 

EMS extended minimum-state 

GAF generalized aerodynamic force 

ISAC interaction structures, aerodynamics, 

and controls 

LS least-squares 

LTI linear time invariant 

MMP modified matrix-Pade 

MS minimum-state 

RFA rational function approximation 
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Table 1. Natural Frequencies of Free-Free 
Symmetric Elastic Modes for ARW-2 


Elastic mode 
number 

Frequency, Hz 

1 

7.85 

2 

14.21 

3 

21.72 

4 

30.27 

5 

33.28 

6 

41.10 

7 

47.01 

8 

63.06 

9 

67.22 

10 

78.24 


Table 2. Baseline (Unoptimized) Approximation Errors 
in Aerodynamic Approximations 


Aerodynamic 

dimension 

Baseline-1, 

J u (a ) 

Baseline-2, 

J c ( b) 

J c /J u 

6 

4.75 

51.88 

10.92 

12 

3.54 

7.36 

2.08 

18 

2.58 

5.90 

2.29 

24 

1.87 

2.59 

1.39 


“Unconstrained approximation error: 


.EE 

' . 2 / ' 
^ Qijd^n) ~~ Qij /Mij 

M 3 i 

n ' 


^Constrained approximation error: 


J c = 


\ 


EE 

3 * 


2 

£ Qij(ik n ) ~ Qijiikn ) / Mij + 2AT. [cj { A m } {j - Cj] 


with constraints given by equations (5.1) (steady state), (5.2) and (5.3) (rigid body 
slopes), and (5.4) (elastic modes near flutter). Capability for explicit inclusion of 
constraints was not utilized. 





Table 3. Optimized, Constrained Approximation Errors in Aerodynamic Approximations 


Aerodynamic 

dimension 

Baseline-2, 

J c (a) 

ELS 

J c 

EMMP 

J c 

EMMP/v 

J c 

EMS 

J c 

4 





9.849 

5 






6 

51.883 

22.57 



4.796 

7 



16.238 

16.238 

4.425 

8 




11.620 

3.442 

9 




9.069 

2.933 

10 




5.760 

2.420 

11 




5.056 


12 

7.362 

4.943 


4.545 


13 




4.151 


14 



4.124 

3.738 


15 




3.302 


16 




2.849 


17 




2.409 


18 

5.900 

3.341 


2.142 


19 




2.015 


20 




1.924 


21 



2.589 

1.833 


22 




1.771 


23 




1.642 


24 

2.591 

2.166 


1.573 


25 




1.517 


26 




1.428 


27 

i 



1.397 


28 



1.828 

1.376 



a Constrained approximation error: 


J c = 


N 


EE 

3 * 


£ 


Qijiikn) Qij{ik n ) 


M ij +2\l[cJ{A r 


’ } ij ^'] 


with constraints given by equations (5.1) (steady state), (5.2) and (5.3) (rigid body slopes), and (5.4) 
(elastic modes near flutter) . Capability for explicit inclusion of constraints was not utilized. 


Table 4. Column Error 



Column 

ELS-24 

EMMP-21 

EMS-10 

EMMP/v- 17 

-plunge 

i 

0.727 

0.688 

0.856 

0.688 

pitch 

2 

.385 

.328 

.415 

.879 

elastic mode 1 

3 

.204 

.248 

.398 

.633 

elastic mode 4 

4 

.504 

.540 

.641 

.540 

elastic mode 6 

5 

.348 

.233 

.399 

.643 

stabilizer 

6 

.970 

.585 

1.300 

.932 

gust 

7 

1.630 

2.318 

1.589 

1.601 
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Qij(i k n ) 


K n 

Real Imaginary 

Unsteady 

0.0 

— 

aerodynamics ► 

0.05 

— 

table 

0.10 

— 


0.20 

— 


0.30 

— 


V 



Figure 1. Approximation of aerodynamic data for oscillatory motion by function defined in complex Laplace 
p-plane. 
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Unsteady 

aerodynamics 

table 


k n 

Qij(ik n ) 

Real Imaginary 

0.0 

-- 

0.05 

— 

0.10 

-- 

0.20 

-- 

0.30 

-- 


Blowup of 
shaded area 



Unconstrained 



Imaginary 



( ik n) 

Q jj (ik) 


Approximation 

error 

eij(ik n ) 


(a) Unconstrained. 

Figure 2. Approximation, with and without constraints, of aerodynamic data for oscillatory motion by function 
defined in complex Laplace p-plane. 
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Unsteady 

aerodynamics 

table 


k n 

Q ij(i k n) 

Real Imaginary 

0.0 

— 

0.05 

-- 

0.10 

-- 

0.20 

— 

0.30 

" 


Blowup of 
shaded area 


n 

Constrained 


V 



Real 


(b) Constrained (at k = 0). 
Figure 2. Concluded. 




Set design 
specifications 

• Order of fit 

• Convergence 

criteria 

• Constraints 



Select initial 

Mo 


Linear 

Minimization 


****** For fixed {b n } *X*!X 

AAA JC * A A A 


Determine all 

^ m ) ij 

which meet 
equality 

constraints 
Evaluate nonlinear 
objective, J c 


Constrained Nonlinear 
Minimization 


Determine new 


Minimum 

tC 

\ i / 


Satisfied 

? 


(a) Extended least-squares (ELS), column- independent formulation. 


Figure 3. Flow of optimization procedures for extended approximation methods. 
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Set design 
specifications 

• Order of fit 

• Convergence 

criteria 

• Constraints 


Start 


Select initial 

(bij } 0 


Linear vC 

** Minimi zation X/* 

A A A 

A ^ A 

AAAAAAAAAAAAAAAA 
AAAAAAAAAAAAAAA 

y For fixed {b^j} y\ 

Determine all 
^ ij 

which meet 
equality 
constraints 
Evaluate nonlinear 
objective, Jj c 


Constrained Nonlinear 
Minimization 


Determine new 


Minimum 

Ji c 


Satisfied 

? 


(b) Extended modified matrix-Pade (EMMP), column-dependent formulation. 


Figure 3. Continued. 
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( 


(c) Extended minimum-state (E 
Figure 3. Conclud 
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Constrained Nonlinear 
Minimization 



formulation. 













(a) Paneling of lifting surfaces for unsteady generalized aerodynamic force calculations. 



Figure 4. Aerodynamic paneling and locations of sensors and control for DAST ARW-2 used in frequency 
response analyses. 
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Total approximation error 



■ Unconstrained 
E2 Constrained 


Aerodynamic dimension 

Figure 6. Unoptimized (baseline) least-squares total approximation errors for unconstrained and constrained 
approximations. 
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Total approximation error 



Figure 7. Optimized, constrained total approximation errors in three extended approaches to approximating 
unsteady aerodynamic forces. 
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Imaginary Imaginary 


1500 



O Tabular data 



(a) Lift due to plunge. 


(b) Pitching moment due to plunge. 




(c) Lift due to pitch. 


(d) Pitching moment due to pitch. 




(e) Lift due to elastic mode 1. (f) Pitching moment due to elastic mode 1. 

Figure 8. RFA plots for selected, optimized ELS, EMMP, and EMS approximations. 



Imaginary Imaginary 






O Tabular data 

Interpolated curve 

♦ ELS-24 

▲ EMMP-2I 

■ EMS-10 



5000 15000 25000 

Real 


(m) Lift due to gust. 



Figure 8. Concluded. 




Tabular data 



(b) Pitching moment due to plunge. 




(c) Lift due to pitch. (d) Pitching moment due to pitch. 




(e) Lift due to elastic mode 1. (f) Pitching moment due to elastic mode 1. 

Figure 9. RFA plots for selected, optimized EMMP approximations for the 17-state case (EMMP/v-17). 
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Imaginary Imaginary Imaginary 








Imaginary Imaginary Imaginary 


O 


Tabular data 
Interpolated curve 
■ EMS/r-10 (relaxed flutter constraints) 








(e) Lift due to elastic mode 1. (f) Pitching moment due to elastic mode 1. 

Figure 10. RFA plots for selected, optimized EMS approximations for the 10-state case in which flutter 
frequency constraints are removed for last two modes. 


47 
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Tabular data 

— Interpolated curve 

- EMS/r-10 (relaxed flutter constraints) 



(g) Lift due to elastic mode 4. 



5 10 15 20 25 

Real 

(h) Pitching moment due to elastic mode 4 



. 


00 -4000 -3500 -3000 -2500 -2000 


Real 



(i) Lift due to elastic mode 6. 


(j) Pitching moment due to elastic mode 6. 
Figure 10. Continued. 





Phase, deg Amplitude, dB 


Interpolated (no rational approximations) 
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Frequency, rad/sec 

Figure 11. Comparison of pitch frequency response due to stabilizer deflection. 
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Figure 12. Comparison of c.g. acceleration frequency response due to stabilizer deflection. 
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